#include <algorithm>
#include <fstream>
#include <strstream>

#include <kuai/misc/strtool.hpp>
#include <kuai/misc/filter.hpp>
#include <kuai/mol/mol.hpp>

namespace kuai {

	std::vector<Element> Element::elements;
	std::map<String, Index> Element::sym2element;

	void Element::update(const Char filename[]) {
		std::ifstream stream(filename);
		update(stream);
	}

	void Element::update() {
		static const Char* szDefaultPeriodicTable = 
		"!#	S(u	OS	VR	IP (ev	EA (ev	A Weight	Cov D	Vdw D	Color	Name	HS	MaxOS\n"
		"0	X	0	0	0	0	0	0	0	0	NONE	0	0	\n"
		"1	H	1	0.32	13.5984	0.7542	1.00794	80	275	1	HYDROGEN	1	1	\n"
		"2	He	0	0.93	24.5874	0	4.0026	400	550	5	HELIUM	0	0	\n"
		"3	Li	1	1.23	5.3917	0.618	6.941	170	305	6	LITHIUM	1	1	\n"
		"4	Be	2	0.9	9.3226	0	9.01218	88	157	5	BERYLLIUM	2	2	\n"
		"5	B	3	0.82	8.298	0.277	10.811	208	387	4	BORON	3	3	\n"
		"6	C	4	0.77	11.2603	1.2619	12.011	180	387	0	CARBON	4	4	\n"
		"7	N	3	0.75	14.5341	1.30	14.0067	170	350	3	NITROGEN	3	5	\n"
		"8	O	2	0.73	13.6181	1.4611	15.9994	170	337	2	OXYGEN	2	2	\n"
		"9	F	1	0.72	17.4228	3.4012	18.9984	160	325	8	FLUORINE	1	1	\n"
		"10	Ne	0	0.71	21.5645	0	20.1797	280	505	5	NEON	0	0	\n"
		"11	Na	1	1.54	5.1391	0.5479	22.9898	243	550	6	SODIUM	1	1	\n"
		"12	Mg	2	1.36	7.6462	0	24.305	275	375	5	MAGNESIUM	2	2	\n"
		"13	Al	3	1.18	5.9858	0.441	26.9815	338	375	9	ALUMINIUM	3	3	\n"
		"14	Si	4	1.11	8.1517	1.385	28.0855	300	550	6	SILICON	4	4	\n"
		"15	P	5	1.06	10.4867	0.7465	30.9738	259	470	6	PHOSPHORUS	3	5	\n"
		"16	S	6	1.02	10.36	2.0771	32.066	255	452	4	SULPHUR	2	6	\n"
		"17	Cl	1	0.99	12.9676	3.6127	35.4527	250	437	8	CHLORINE	1	7	\n"
		"18	Ar	0	0.98	15.7596	0	39.948	392	692	5	ARGON	0	0	\n"
		"19	K	1	2.03	4.3407	0.5015	39.0983	332	597	6	POTASSIUM	1	1	\n"
		"20	Ca	2	1.74	6.1132	0.0184	40.078	248	487	5	CALCIUM	2	2	\n"
		"21	Sc	3	1.44	6.5614	0.188	44.9559	360	330	9	SCANDIUM	0	3	\n"
		"22	Ti	4	1.32	6.8282	0.079	47.88	368	487	9	TITANIUM	0	4	\n"
		"23	V	5	1.22	6.7463	0.525	50.9415	332	265	9	VANADIUM	0	5	\n"
		"24	Cr	6	1.18	6.7666	0.666	51.996	338	282	9	CHROMIUM	0	6	\n"
		"25	Mn	4	1.17	7.434	0	54.938	338	297	9	MANGANESE	0	6	\n"
		"26	Fe	3	1.17	7.9024	0.151	55.847	335	487	8	IRON	0	6	\n"
		"27	Co	2	1.16	7.881	0.662	58.9332	332	282	9	COBALT	0	3	\n"
		"28	Ni	2	1.15	7.6398	1.156	58.6934	405	310	6	NICKEL	0	3	\n"
		"29	Cu	2	1.17	7.7264	1.235	63.546	380	287	6	COPPER	0	2	\n"
		"30	Zn	2	1.25	9.3941	0	65.39	362	287	6	ZINC	0	2	\n"
		"31	Ga	3	1.26	5.9993	0.3	69.723	305	387	9	GALLIUM	3	3	\n"
		"32	Ge	4	1.22	7.9	1.233	72.61	292	999	9	GERMANIUM	4	4	\n"
		"33	As	3	1.2	9.8152	0.81	74.9216	302	207	9	ARSENIC	3	5	\n"
		"34	Se	4	1.16	9.7524	2.0207	78.96	305	225	9	SELENIUM	2	6	\n"
		"35	Br	1	1.14	11.8138	3.3636	79.904	302	437	8	BROMINE	1	7	\n"
		"36	Kr	0	1.89	13.9996	0	83.8	400	475	9	KRYPTON	0	0	\n"
		"37	Rb	1	2.16	4.1771	0.4859	85.4678	368	662	9	RUBIDIUM	1	1	\n"
		"38	Sr	2	1.91	5.6948	0.11	87.62	280	505	5	STRONTIUM	2	2	\n"
		"39	Y	3	1.62	6.217	0.307	88.9059	445	402	9	YTTRIUM	0	3	\n"
		"40	Zr	4	1.45	6.6339	0.426	91.224	390	355	9	ZIRCONIUM	0	4	\n"
		"41	Nb	5	1.34	6.7589	0.893	92.9064	370	332	9	NIOBIUM	0	5	\n"
		"42	Mo	6	1.3	7.0924	0.746	95.94	368	437	9	MOLYBDENUM	0	6	\n"
		"43	Tc	7	1.27	7.28	0.55	98	338	450	9	TECHNETIUM	0	7	\n"
		"44	Ru	4	1.25	7.3605	1.05	101.07	350	300	9	RUTHENIUM	0	6	\n"
		"45	Rh	3	1.25	7.4589	1.137	102.906	362	305	9	RHODIUM	0	4	\n"
		"46	Pd	2	1.28	8.3369	0.557	106.42	375	360	9	PALLADIUM	0	4	\n"
		"47	Ag	1	1.34	7.5762	1.302	107.868	398	387	9	SILVER	0	1	\n"
		"48	Cd	2	1.41	8.9937	0	112.41	422	437	9	CADMIUM	0	2	\n"
		"49	In	3	1.44	5.7864	0.3	114.82	408	362	9	INDIUM	3	3	\n"
		"50	Sn	4	1.41	7.3438	0.112	118.71	365	417	9	TIN	4	4	\n"
		"51	Sb	3	1.4	8.64	1.07	121.757	365	280	9	ANTIMONY	3	5	\n"
		"52	Te	4	1.36	9.0096	1.9708	127.6	368	315	9	TELLURIUM	2	6	\n"
		"53	I	1	1.33	10.4513	3.059	126.905	350	437	8	IODINE	1	7	\n"
		"54	Xe	0	1.31	12.1299	0	131.29	425	525	9	XENON	0	0	\n"
		"55	Cs	1	2.35	3.8939	0.4716	132.905	418	752	9	CAESIUM	1	1	\n"
		"56	Ba	2	1.98	5.2117	0.15	137.33	335	602	8	BARIUM	2	2	\n"
		"57	La	3	1.25	5.577	0.5	138.905	468	457	9	LANTHANUM	0	3	\n"
		"58	Ce	4	1.65	5.5387	0	140.12	458	465	9	CERIUM	0	4	\n"
		"59	Pr	4	1.65	5.464	0	140.908	455	405	9	PRASEODYMIUM	0	4	\n"
		"60	Nd	3	1.64	5.525	0	144.24	452	447	9	NEODYMIUM	0	3	\n"
		"61	Pm	3	1.63	5.55	0	145	450	440	9	PROMETHIUM	0	3	\n"
		"62	Sm	3	1.62	5.6437	0	150.36	450	435	9	SAMARIUM	0	3	\n"
		"63	Eu	3	1.85	5.6704	0	151.965	498	490	9	EUROPIUM	0	3	\n"
		"64	Gd	3	1.61	6.15	0	157.25	448	422	9	GADOLINIUM	0	3	\n"
		"65	Tb	3	1.59	5.8639	0	158.925	440	415	9	TERBIUM	0	4	\n"
		"66	Dy	3	1.59	5.9389	0	162.5	438	407	9	DYSPROSIUM	0	3	\n"
		"67	Ho	3	1.58	6.0216	0	164.93	435	402	9	HOLMIUM	0	3	\n"
		"68	Er	3	1.57	6.1078	0	167.26	432	397	9	ERBIUM	0	3	\n"
		"69	Tm	3	1.56	6.1843	0	168.934	430	392	9	THULIUM	0	3	\n"
		"70	Yb	3	1.7	6.2542	0	173.04	485	385	9	YTTERBIUM	0	3	\n"
		"71	Lu	3	1.56	5.4259	0	174.967	430	382	9	LUTETIUM	0	3	\n"
		"72	Hf	4	1.44	6.8251	0	178.49	392	350	9	HAFNIUM	0	4	\n"
		"73	Ta	5	1.34	7.89	0.322	180.948	358	305	9	TANTALUM	0	5	\n"
		"74	W	6	1.3	7.98	0.815	183.85	342	315	9	TUNGSTEN	0	6	\n"
		"75	Re	7	1.28	7.88	0.15	186.207	338	325	9	RHENIUM	0	7	\n"
		"76	Os	4	1.26	8.7	1.1	190.2	342	395	9	OSMIUM	0	6	\n"
		"77	Ir	4	1.27	9.1	1.565	192.22	330	305	9	IRIDIUM	0	6	\n"
		"78	Pt	4	1.3	9	2.128	195.08	375	387	9	PLATINUM	0	4	\n"
		"79	Au	3	1.34	9.2257	2.3086	196.967	375	362	6	GOLD	0	3	\n"
		"80	Hg	2	1.49	10.4375	0	200.59	425	495	9	MERCURY	0	2	\n"
		"81	Tl	1	1.48	6.1083	0.2	204.282	388	427	9	THALLIUM	3	3	\n"
		"82	Pb	2	1.47	7.4167	0.364	207.2	385	540	9	LEAD	4	4	\n"
		"83	Bi	3	1.46	7.289	0.946	208.98	385	432	9	BISMUTH	3	5	\n"
		"84	Po	4	1.53	8.4167	1.9	209	420	302	9	POLONIUM	2	4	\n"
		"85	At	1	1.47	10	2.8	210	302	280	9	ASTATINE	1	7	\n"
		"86	Rn	0	1.5	10.7485	0	222	475	575	9	RADON	0	0	\n"
		"87	Fr	1	2.4	3.5	0.46	223	450	810	9	FRANCIUM	1	1	\n"
		"88	Ra	2	2	5.2789	0	226.025	358	642	9	RADIUM	2	2	\n"
		"89	Ac	3	1.5	5.17	0	227	295	530	9	ACTINIUM	0	3	\n"
		"90	Th	4	1.65	6.08	0	232.038	255	460	9	THORIUM	0	4	\n"
		"91	Pa	5	1.65	5.89	0	231.036	222	400	9	PROTACTINIUM	0	5	\n"
		"92	U	6	1.42	6.1941	0	238.029	242	437	9	URANIUM	0	6	\n"
		"93	Np	6	1.4	6.2657	0	237.048	238	427	9	NEPTUNIUM	0	6	\n"
		"94	Pu	6	1.08	6.06	0	244	232	417	9	PLUTONIUM	0	6	\n"
		"95	Am	6	1.5	5.993	0	243	230	415	9	AMERICIUM	0	6	\n"
		"96	Cm	3	1.5	6.02	0	247	228	412	9	CURIUM	0	3	\n"
		"97	Bk	4	1.5	6.23	0	247	225	410	9	BERKELIUM	0	4	\n"
		"98	Cf	3	1.5	6.3	0	251	222	407	9	CALIFORNIUM	0	3	\n"
		"99	Es	3	1.5	6.42	0	252	220	405	9	EINSTEINIUM	0	3	\n"
		"100	Fm	3	1.5	6.5	0	257	218	402	9	FERMIUM	0	3	\n"
		"101	Md	3	1.5	6.58	0	258	215	400	9	MENDELEVIUM	0	3	\n"
		"102	No	3	1.5	6.65	0	259	212	397	9	NOBELIUM	0	1	\n"
		"103	Lr	3	1.5	6.6	0	260	210	395	9	LAWRENCIUM	0	1	\n"
		"";

		std::istrstream stream(szDefaultPeriodicTable);
		update(stream);
	}

	namespace {
		bool readElement(const String& line, Element& v0)
		{
			StringArray tokens = split(line);
			try { 
				v0.number = lexical_cast<size_t>(tokens[0]);
				v0.symbol = tokens[1];
				v0.commonValence = lexical_cast<size_t>(tokens[2]);
				v0.valenceRadii = lexical_cast<RealNumber>(tokens[3]);
				v0.ionPotential = lexical_cast<RealNumber>(tokens[4]);
				v0.elecAffinities = lexical_cast<RealNumber>(tokens[5]);
				v0.weight = lexical_cast<RealNumber>(tokens[6]);
				// v0.cov_d = boost::lexical_cast<int>(tokens[7]);
				// v0.vdw_d = boost::lexical_cast<int>(tokens[8]);
				// v0.color = boost::lexical_cast<size_t>(tokens[9]);
				// v0.name = tokens[10];
				v0.implicitH = lexical_cast<size_t>(tokens[11]);
				v0.maxValence = lexical_cast<size_t>(tokens[12]);
				return true;
			}
			catch(boost::bad_lexical_cast&)
			{ }
			return false;
		}
	}

	void Element::update(std::istream& stream) {
		String line;
		Element e;
		elements.reserve(128);
		while (getline(stream, line)) {
			boost::trim(line);
			if (isdigit(line[0])) {
				if (readElement(line, e)) {
					if (elements.size() <= e.number) {
						elements.resize(e.number+1);
					}
					elements[e.number] = e;
					sym2element[e.symbol] = e.number;
				}
			}
		}
	}

	const Element& Element::getElement(Index i) {
		if (elements.empty()) {
			update();
		}
		if (i >= elements.size()) {
			i = 0;
		}
		return elements[i];
	}

	const Element& Element::getElement(const Char symbol[]) {
		return getElement(String(symbol));
	}

	const Element& Element::getElement(const String& symbol) {
		if (elements.empty()) {
			update();
		}
		Index i = 0;
		std::map<String, Index>::const_iterator it = sym2element.find(symbol);
		if (it != sym2element.end()) {
			i = it->second;
		}
		return getElement(i);
	}

	extern const Char BOND_SYMBOLS[8] = " .-:=*#";

	void guessBondLink(Index nAtoms, const Index atoms[], const XYZ coords[], const PBC* pbc, 
		std::vector<Index>& bonds, std::map<Index, RealNumber>& valenceRadii, RealNumber tor)
	{
		for (size_t i = 0; i < nAtoms; ++i) {
			if (valenceRadii.find(atoms[i]) == valenceRadii.end()) {
				valenceRadii.insert(std::make_pair(atoms[i], Element::getElement(atoms[i]).valenceRadii));
			}
		}
		LowerTriangleFilter0 filter(nAtoms, true);
		for (size_t i = 0; i < bonds.size(); i += 3)
		{
			filter.reset(bonds[i+1], bonds[i+2]);
		}
		bonds.reserve(nAtoms*3);
		for (size_t i = 0; i < nAtoms; ++i) {
			RealNumber rI = valenceRadii[atoms[i]];
			const XYZ& coordI=coords[i];
			for (size_t j = i+1; j < nAtoms; ++j) {
				if (filter.test(i, j)) {
					const XYZ& coordJ=coords[j];
					XYZ delta(coordI.x-coordJ.x, coordI.y-coordJ.y, coordI.z-coordJ.z);
					if (pbc) {
						pbc->update(delta);
					}
					RealNumber rr = rI + valenceRadii[atoms[j]] + tor;
					if (shorter(delta, rr)) {
						bonds.push_back(Bond::SINGLE_BOND);
						bonds.push_back(i);
						bonds.push_back(j);
					}
				}
			}
		}
	}


	RealNumber formulaWeight(const String& formula) {
		RealNumber result=0;
		
		for (const char* p = formula.c_str(); *p != NULL; ) {
			size_t factor=1;
			RealNumber part = 0;
			const char* q;
			for (q = p; isdigit(*q); ++q) 
			{ }
			if (p != q) {
				factor = lexical_cast<size_t>(String(p, q));
			}
			p = q;
			while (*p != '.' && *p != char(0)) {
				q = p + 1;
				if (islower(*q)) {
					q = q + 1;
				}
				String symbol = String(p, q);
				size_t n = 1;
				p = q;
				while (isdigit(*q)) {
					++q;
				}
				if (p != q) {
					n = lexical_cast<size_t>(String(p, q));
					p = q;
				}
				part += Element::getElement(symbol.c_str()).weight * n;
			}
			result += part * factor;
			while (*p == '.') {
				++p;
			}
		}
		return result;
	}

}
